function R = median_features(image)
if size(image,3) == 3
    image = rgb2gray(image);
end
f3_1 = processOne(image, [1 0 1;0 1 0;1 0 1]);
f3_2 = processOne(image, [0 1 0;1 1 1;0 1 0]);
f5_1 = processOne(image, [0 1 0 1 0;1 1 0 1 1;0 0 1 0 0;1 1 0 1 1;0 1 0 1 0]);
f5_2 = processOne(image, [1 0 1 0 1;0 0 0 0 0;1 0 1 0 1;0 0 0 0 0;1 0 1 0 1]);
f5_3 = processOne(image, [0 0 1 0 0;0 1 1 1 0;1 1 1 1 1;0 1 1 1 0;0 0 1 0 0]);
f7_1 = processOne(image, [0 1 1 0 1 1 0;1 1 1 0 1 1 1;1 1 1 0 1 1 1;0 0 0 1 0 0 0;1 1 1 0 1 1 1;1 1 1 0 1 1 1;0 1 1 0 1 1 0]);
f7_2 = processOne(image, [0 0 0 1 0 0 0;0 1 0 1 0 1 0;0 0 0 1 0 0 0;1 1 1 1 1 1 1;0 0 0 1 0 0 0;0 1 0 1 0 1 0;0 0 0 1 0 0 0]);
f7_3 = processOne(image, [1 0 1 0 1 0 1;0 0 0 0 0 0 0;1 0 1 0 1 0 1;0 0 0 1 0 0 0;1 0 1 0 1 0 1;0 0 0 0 0 0 0;1 0 1 0 1 0 1]);
f7_4 = processOne(image, [0 0 0 1 0 0 0;0 0 1 1 1 0 0;0 1 1 1 1 1 0;1 1 1 1 1 1 1;0 1 1 1 1 1 0;0 0 1 1 1 0 0;0 0 0 1 0 0 0]);
R = [f3_1 f3_2 f5_1 f5_2 f5_3 f7_1 f7_2 f7_3 f7_4]; 
end

function f = processOne(image, filter)
N = (sum(filter(:)) + 1)/2;
r = ordfilt2(image, N, filter) - image;
%MAPPING = getmapping(8, 'u2');
f = [lpq(r, 3, 0, 1, 'nh'), ...
    lpq(r, 5, 0, 1, 'nh'), ...
    lpq(r, 7, 0, 1, 'nh')];
end
